#include <iostream>
#include <assert.h>
#include "MeshInfoProcess2D.hpp"


MeshInfo_2d::MeshInfo_2d( vector<double> sKnots,
    vector<double> tKnots,
    int s_Degree,
    int t_Degree,
    int nqpx,
    int nqpy)
{
  sDegree = s_Degree;
  tDegree = t_Degree;
  nFunc_x = sKnots.size() - sDegree - 1;
  nFunc_y = tKnots.size() - tDegree - 1;
  nFunc   = nFunc_x * nFunc_y;
  nElem_x = nFunc_x - sDegree;
  nElem_y = nFunc_y - tDegree;
  nElem   = nElem_x * nElem_y;
  nLocBas = (sDegree+1) * (tDegree+1);
  nqp_x   = nqpx;
  nqp_y   = nqpy;

  int ii;
  hx = new double [nElem_x];
  hy = new double [nElem_y];

  for(ii=0; ii<nElem_x; ii++)
    hx[ii] = sKnots[sDegree+ii+1] - sKnots[sDegree+ii];

  for(ii=0; ii<nElem_y; ii++)
    hy[ii] = tKnots[tDegree+ii+1] - tKnots[tDegree+ii];


  int ** IENv;
  IEN  = new int * [nElem];
  IENv = new int * [nElem];
  INN  = new int * [nFunc];

  for(ii=0; ii<nElem; ii++)
  {
    IENv[ii] = new int [nLocBas];
    IEN[ii] = new int [nLocBas];
  }
  for(ii=0; ii<nFunc; ii++)
    INN[ii] = new int [2];

  int A = -1;
  int e = -1;
  int B = 0;
  int b = 0;
  for(int j=0; j<nFunc_y; j++)
  {
    for(int i=0; i<nFunc_x; i++)
    {
      A = A + 1;
      INN[A][0] = i;
      INN[A][1] = j;
      if( (i>=sDegree) && (j>=tDegree) )
      {
        e = e+1;
        for(int jloc = 0; jloc<=tDegree; jloc++)
        {
          for(int iloc=0; iloc<=sDegree; iloc++)
          {
            B = A - jloc * nFunc_x - iloc;
            b = jloc * (sDegree + 1) + iloc;
            IENv[e][b] = B;
          }
        }
      }
    }
  }

  for(int i=0; i<nElem; i++)
  {
    for(int j=0; j<nLocBas; j++)
      IEN[i][j] = IENv[i][nLocBas-1-j];
  }

  for(int i=0; i<nElem; i++)
    delete [] IENv[i];
      delete [] IENv;

  qp_x = new double[nqp_x];
  wq_x = new double[nqp_x];
  qp_y = new double[nqp_y];
  wq_y = new double[nqp_y];

  QuadratureRule(nqp_x, qp_x, wq_x);
  QuadratureRule(nqp_y, qp_y, wq_y);
};

MeshInfo_2d::~MeshInfo_2d( )
{
  for(int i=0; i<nElem; i++)
    delete [] IEN[i];
  for(int i =0; i<nFunc; i++)
    delete [] INN[i];
  delete [] IEN;
  delete [] INN;
  delete [] hx;
  delete [] hy;
  delete [] qp_x;
  delete [] qp_y;
  delete [] wq_x;
  delete [] wq_y;
  std::cout<<"CLEAN: MeshInfo2d class deleted."<<std::endl;
};

void MeshInfo_2d::QuadratureRule(int nqp, double * &qp, double * &wq)
{
  switch (nqp)
  {
    case 1:
      qp[0] = 0.5;
      wq[0] = 1.0;
      break;
    case 2:
      qp[0] = 0.788675134594813;
      wq[0] = 0.5;
      qp[1] = 0.211324865405187;
      wq[1] = 0.5;
      break;
    case 3:
      qp[0] = 0.887298334620742;
      qp[1] = 0.5;
      qp[2] = 0.112701665379258;
      wq[0] = 0.277777777777777;
      wq[1] = 0.444444444444444;
      wq[2] = 0.277777777777777;
      break;
    case 4:
      qp[0] = 0.930568155797026;
      qp[1] = 0.669990521792428;
      qp[2] = 0.330009478207572;
      qp[3] = 0.069431844202974;

      wq[0] = 0.173927422568727;
      wq[1] = 0.326072577431273;
      wq[2] = 0.326072577431273;
      wq[3] = 0.173927422568727;
      break;
    default:
      cerr<<"ERROR: We have quadrature rule with up to 4 quadrature points"<<endl;
      assert(0);
  }
}

